library(readxl)
library(ppcor) # Para correlaciones parciales
library(dplyr)
library(ggplot2)
library(fastDummies) # Para crear variables dummy7 Regresion lineal multiple
7.1 Introducción
Este documento presenta una guía detallada para el análisis de Regresión Lineal Múltiple a través de dos ejercicios prácticos. El primer ejercicio se centra en un modelo con predictores continuos, mientras que el segundo explora un modelo con predictores categóricos (ANOVA/ANCOVA). Se abordan todos los pasos del modelado: análisis descriptivo, estimación de parámetros, pruebas de hipótesis y construcción de intervalos de confianza.
7.1.1 Librerías Requeridas
8 EJERCICIO 1: Masa Corporal Magra en Atletas
Se busca modelar la masa corporal magra (lbm) en función de la altura (ht), el peso (wt) y el recuento de glóbulos rojos (rcc).
8.1 a) Análisis Descriptivo
Se inicia con una exploración de las relaciones entre las variables mediante diagramas de dispersión y matrices de correlación.
Atletas <- read_excel("Atletas.xlsx")
# View(Atletas)
# Matriz de diagramas de dispersión
plot(Atletas)
# Matriz de correlación de Pearson
cor(Atletas) lbm ht wt rcc
lbm 1.00000000 0.71132706 0.93917965 0.08524203
ht 0.71132706 1.00000000 0.71506428 0.01460278
wt 0.93917965 0.71506428 1.00000000 0.02054923
rcc 0.08524203 0.01460278 0.02054923 1.00000000# Matriz de correlaciones parciales
pcor(Atletas)$estimate
lbm ht wt rcc
lbm 1.0000000 0.1687532 0.8798006 0.1947647
ht 0.1687532 1.0000000 0.1865226 -0.0329934
wt 0.8798006 0.1865226 1.0000000 -0.1646123
rcc 0.1947647 -0.0329934 -0.1646123 1.0000000
$p.value
lbm ht wt rcc
lbm 0.000000e+00 0.1002613 4.047452e-32 0.05722977
ht 1.002613e-01 0.0000000 6.881920e-02 0.74963535
wt 4.047452e-32 0.0688192 0.000000e+00 0.10900297
rcc 5.722977e-02 0.7496354 1.090030e-01 0.00000000
$statistic
lbm ht wt rcc
lbm 0.000000 1.6599295 17.944907 1.9251815
ht 1.659929 0.0000000 1.840707 -0.3200571
wt 17.944907 1.8407070 0.000000 -1.6180485
rcc 1.925182 -0.3200571 -1.618048 0.0000000
$n
[1] 98
$gp
[1] 2
$method
[1] "pearson"8.2 c) Estimación de los Parámetros del Modelo
Se ajusta un modelo de Regresión Lineal Múltiple de la forma: \[ \text{lbm}_i = \beta_0 + \beta_1 \text{ht}_i + \beta_2 \text{wt}_i + \beta_3 \text{rcc}_i + \epsilon_i \] Los parámetros se estiman por Mínimos Cuadrados Ordinarios (MCO).
fit <- lm(lbm ~ 1 + ht + wt + rcc, data = Atletas)
summary(fit)
Call:
lm(formula = lbm ~ 1 + ht + wt + rcc, data = Atletas)
Residuals:
Min 1Q Median 3Q Max
-5.0163 -1.8628 0.1932 1.8690 6.0228
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.30797 6.71133 -0.195 0.8459
ht 0.06848 0.04126 1.660 0.1003
wt 0.56637 0.03156 17.945 <2e-16 ***
rcc 1.42480 0.74008 1.925 0.0572 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.359 on 94 degrees of freedom
Multiple R-squared: 0.8896, Adjusted R-squared: 0.8861
F-statistic: 252.6 on 3 and 94 DF, p-value: < 2.2e-16# Intervalos de confianza para los coeficientes
confint(fit) 2.5 % 97.5 %
(Intercept) -14.63346824 12.0175228
ht -0.01343245 0.1503940
wt 0.50370761 0.6290412
rcc -0.04465812 2.89425168.3 d) Comprobación Matricial de los Parámetros
Los estimadores MCO pueden calcularse manualmente usando la notación matricial: \[ \hat{\boldsymbol{\beta}} = (X^T X)^{-1} X^T Y \]
n <- dim(Atletas)[1]
Y <- as.matrix(Atletas[, 1])
X <- as.matrix(Atletas[, 2:4])
X <- cbind(rep(1, n), X) # Añadir columna de unos para el intercepto
colnames(X)[1]<-"interc"
betahat <- solve(t(X) %*% X) %*% t(X) %*% Y
# Comprobar si los cálculos manuales y los de lm() coinciden
all.equal(as.numeric(betahat), as.numeric(fit$coefficients))[1] TRUE8.4 e) Matriz de Varianzas y Covarianzas Estimada
La matriz de varianza-covarianza de los estimadores se calcula como: \[ \text{Var}(\hat{\boldsymbol{\beta}}) = \hat{\sigma}^2 (X^T X)^{-1} \] donde \(\hat{\sigma}^2 = \frac{SCE}{n-p}\) es la varianza residual estimada.
# Resultado de R
vcov(fit) (Intercept) ht wt rcc
(Intercept) 45.04189942 -2.346794e-01 0.0971874453 -2.391439e+00
ht -0.23467939 1.701993e-03 -0.0009309838 3.986147e-06
wt 0.09718745 -9.309838e-04 0.0009961501 -3.377628e-04
rcc -2.39143869 3.986147e-06 -0.0003377628 5.477249e-01# Cálculo manual
varhat <- (1 / (n - 4)) * sum((Y - X %*% betahat)^2)
vcovbetahat <- varhat * solve(t(X) %*% X)
# Comprobar si son iguales
all.equal(as.numeric(vcovbetahat), as.numeric(vcov(fit)))[1] TRUE8.5 i) Prueba de Significancia Global del Modelo
Se evalúa si el modelo en su conjunto es significativo. - \(H_0\): \(\beta_1 = \beta_2 = \beta_3 = 0\) (El modelo no es significativo). - \(H_1\): Al menos un \(\beta_j \neq 0\).
# Opción 1: Usar el F-statistic de la salida de summary()
summary(fit)
Call:
lm(formula = lbm ~ 1 + ht + wt + rcc, data = Atletas)
Residuals:
Min 1Q Median 3Q Max
-5.0163 -1.8628 0.1932 1.8690 6.0228
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.30797 6.71133 -0.195 0.8459
ht 0.06848 0.04126 1.660 0.1003
wt 0.56637 0.03156 17.945 <2e-16 ***
rcc 1.42480 0.74008 1.925 0.0572 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.359 on 94 degrees of freedom
Multiple R-squared: 0.8896, Adjusted R-squared: 0.8861
F-statistic: 252.6 on 3 and 94 DF, p-value: < 2.2e-16# Opción 2: Comparar el modelo completo con un modelo reducido (solo intercepto) usando anova()
fit0 <- lm(lbm ~ 1, data = Atletas)
anova(fit0, fit, test = "F")Analysis of Variance Table
Model 1: lbm ~ 1
Model 2: lbm ~ 1 + ht + wt + rcc
Res.Df RSS Df Sum of Sq F Pr(>F)
1 97 4741.4
2 94 523.2 3 4218.2 252.59 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusión: En ambas opciones, el p-valor es extremadamente pequeño, por lo que se rechaza \(H_0\). El modelo es globalmente significativo.
8.6 m) Prueba de Significancia de un Subconjunto de Variables
Se prueba si las variables ht y rcc son conjuntamente necesarias en el modelo. - \(H_0\): \(\beta_{ht} = \beta_{rcc} = 0\). - \(H_1\): Al menos uno de los dos coeficientes es distinto de cero.
# Se compara el modelo completo con un modelo reducido que excluye ht y rcc
fit1 <- lm(lbm ~ 1 + wt, data = Atletas)
anova(fit1, fit, test = "F")Analysis of Variance Table
Model 1: lbm ~ 1 + wt
Model 2: lbm ~ 1 + ht + wt + rcc
Res.Df RSS Df Sum of Sq F Pr(>F)
1 96 559.21
2 94 523.24 2 35.964 3.2304 0.04397 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusión: El p-valor pequeño indica que se rechaza \(H_0\). Las variables ht y rcc son conjuntamente significativas y deben permanecer en el modelo.
8.7 n) Prueba de Hipótesis Lineal: \(\beta_{ht} = \beta_{wt}\) ?
Se prueba si los efectos de la altura y el peso son iguales. - \(H_0\): \(\beta_{ht} - \beta_{wt} = 0\). - \(H_1\): \(\beta_{ht} - \beta_{wt} \neq 0\).
a <- c(0, 1, -1, 0) # Vector para la combinación lineal
t_C <- t(a)%*%as.matrix(fit$coeff)/sqrt((t(a)%*%as.matrix(vcov(fit))%*%a))
pval_n <- 2 * (1 - pt(abs(t_C), n - 4))
c(t_calculado = t_C, p_valor = pval_n) t_calculado p_valor
-7.373079e+00 6.449596e-11 Conclusión: Se rechaza \(H_0\). Los efectos de la altura y el peso son significativamente diferentes.
8.8 o) Intervalo de Confianza del 95% para la Media
Se calcula un intervalo de confianza para el valor esperado de lbm para un atleta con ht=170, wt=70 y rcc=4.5.
nuevos_datos <- data.frame(ht = 170, wt = 70, rcc = 4.5)
predict(fit, nuevos_datos, se.fit = TRUE, interval = "confidence", level = 0.95)$fit
fit lwr upr
1 56.39155 55.67493 57.10817
$se.fit
[1] 0.3609235
$df
[1] 94
$residual.scale
[1] 2.3593279 EJERCICIO 2: Índice de Placa Bacteriana (IPB)
Se busca modelar el IPB en función del grupo de tratamiento y el sexo del paciente, un ejemplo de modelo ANCOVA.
9.1 a) Estadísticas Descriptivas
ipb_data <- read_excel("ipb.xlsx")
head(ipb_data)# A tibble: 6 × 3
grupo sexo ipb
<chr> <chr> <dbl>
1 Grupo1 F 1.2
2 Grupo1 F 1.43
3 Grupo1 F 1.1
4 Grupo1 F 1.45
5 Grupo1 F 0.95
6 Grupo1 F 2.75# Resúmenes por grupo
ipb_data %>%
group_by(grupo, sexo) %>%
summarise(
n = n(),
media_ipb = mean(ipb),
sd_ipb = sd(ipb)
)`summarise()` has grouped output by 'grupo'. You can override using the
`.groups` argument.# A tibble: 6 × 5
# Groups: grupo [3]
grupo sexo n media_ipb sd_ipb
<chr> <chr> <int> <dbl> <dbl>
1 Grupo1 F 10 1.35 0.541
2 Grupo1 M 10 2.72 0.527
3 Grupo2 F 10 1.51 0.570
4 Grupo2 M 10 3.45 0.476
5 Grupo3 F 10 1.89 0.624
6 Grupo3 M 10 3.82 0.420# Boxplot
boxplot(ipb ~ grupo * sexo, data = ipb_data,
xlab = "Grupo y Sexo", ylab = "IPB",
frame = FALSE, col = c("#00AFBB", "#E7B800", "#FC4E07"))
9.2 c) Estimación de Parámetros del Modelo
Se ajusta un modelo con predictores categóricos. R crea automáticamente las variables dummy. \[ \text{ipb}_i = \beta_0 + \beta_1 \text{grupo2}_i + \beta_2 \text{grupo3}_i + \beta_3 \text{sexoM}_i + \epsilon_i \] - \(\beta_0\) es la media del grupo de referencia (Grupo 1, Femenino). - Los otros \(\beta\) son los cambios promedio con respecto al grupo de referencia.
fit_ipb <- lm(ipb ~ 1 + grupo + sexo, data = ipb_data)
summary(fit_ipb)
Call:
lm(formula = ipb ~ 1 + grupo + sexo, data = ipb_data)
Residuals:
Min 1Q Median 3Q Max
-1.22478 -0.29049 0.02787 0.27287 1.58787
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.1621 0.1390 8.361 1.98e-11 ***
grupoGrupo2 0.4464 0.1702 2.622 0.0112 *
grupoGrupo3 0.8225 0.1702 4.832 1.09e-05 ***
sexoM 1.7457 0.1390 12.560 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5383 on 56 degrees of freedom
Multiple R-squared: 0.7639, Adjusted R-squared: 0.7512
F-statistic: 60.38 on 3 and 56 DF, p-value: < 2.2e-169.3 d) ¿El IPB depende del Grupo?
Se prueba si los coeficientes asociados a los grupos son conjuntamente cero. - \(H_0\): \(\beta_{grupo2} = \beta_{grupo3} = 0\).
# Se compara el modelo completo con un modelo reducido sin la variable 'grupo'
fit_reducido_grupo <- lm(ipb ~ 1 + sexo, data = ipb_data)
anova(fit_reducido_grupo, fit_ipb, test = "F")Analysis of Variance Table
Model 1: ipb ~ 1 + sexo
Model 2: ipb ~ 1 + grupo + sexo
Res.Df RSS Df Sum of Sq F Pr(>F)
1 58 23.011
2 56 16.229 2 6.782 11.701 5.674e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusión: El p-valor es muy pequeño, se rechaza \(H_0\). El grupo tiene un efecto significativo en el IPB.
9.4 e) ¿El IPB depende del Sexo?
Se prueba si el coeficiente asociado al sexo es cero. - \(H_0\): \(\beta_{sexoM} = 0\).
# Se puede juzgar directamente por el p-valor de 'sexoM' en summary(fit_ipb)
summary(fit_ipb)
Call:
lm(formula = ipb ~ 1 + grupo + sexo, data = ipb_data)
Residuals:
Min 1Q Median 3Q Max
-1.22478 -0.29049 0.02787 0.27287 1.58787
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.1621 0.1390 8.361 1.98e-11 ***
grupoGrupo2 0.4464 0.1702 2.622 0.0112 *
grupoGrupo3 0.8225 0.1702 4.832 1.09e-05 ***
sexoM 1.7457 0.1390 12.560 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5383 on 56 degrees of freedom
Multiple R-squared: 0.7639, Adjusted R-squared: 0.7512
F-statistic: 60.38 on 3 and 56 DF, p-value: < 2.2e-16Conclusión: El p-valor (0.13) es mayor que 0.05. No hay evidencia de que el sexo tenga un efecto significativo en el IPB, después de controlar por el grupo.
9.5 f) ¿Hay Interacción entre Sexo y Grupo?
Se prueba si el efecto del grupo sobre el IPB es el mismo para hombres y mujeres. - \(H_0\): No hay interacción.
# Se ajusta un modelo con el término de interacción
fit_interaccion <- lm(ipb ~ 1 + grupo + sexo + grupo:sexo, data = ipb_data)
summary(fit_interaccion)
Call:
lm(formula = ipb ~ 1 + grupo + sexo + grupo:sexo, data = ipb_data)
Residuals:
Min 1Q Median 3Q Max
-1.13275 -0.31603 -0.04393 0.24413 1.40200
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.3480 0.1677 8.037 8.61e-11 ***
grupoGrupo2 0.1667 0.2372 0.703 0.4852
grupoGrupo3 0.5446 0.2372 2.296 0.0256 *
sexoM 1.3740 0.2372 5.792 3.66e-07 ***
grupoGrupo2:sexoM 0.5594 0.3355 1.668 0.1012
grupoGrupo3:sexoM 0.5558 0.3355 1.657 0.1034
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5304 on 54 degrees of freedom
Multiple R-squared: 0.7789, Adjusted R-squared: 0.7585
F-statistic: 38.05 on 5 and 54 DF, p-value: < 2.2e-16# Se comparan los modelos con y sin interacción
anova(fit_ipb, fit_interaccion, test = "F")Analysis of Variance Table
Model 1: ipb ~ 1 + grupo + sexo
Model 2: ipb ~ 1 + grupo + sexo + grupo:sexo
Res.Df RSS Df Sum of Sq F Pr(>F)
1 56 16.229
2 54 15.192 2 1.0364 1.8419 0.1683Conclusión: El p-valor (0.55) es grande. No hay evidencia de un efecto de interacción significativo.